/* Copyright (c) Colorado School of Mines, 2011.*/
/* All rights reserved.                       */

/* SUPHASE: $Revision: 1.2 $ ; $Date: 2011/11/20 18:24:54 $    */

#include "su.h"
#include "segy.h"

/******************************** self documentation ******************/
char *sdoc[] = {
   "                                                                  ",
   " SUPHASE - PHASE manipulation                                     ",
   "                                                                  ",
   "  suphase  <stdin >sdout                                          ",
   "                                                                  ",
   " Required parameters:                                             ",
   " none                                                             ",
   " Optional parameters:                                             ",
   " a=90         constant phase shift (degrees)                      ",
   " b=1.0        multiplicative phase shift (???)                    ",
   " c=0.0        time shift (seconds)                                ",
   "                                                                  ",
   "phase = a + b*(old_phase) + c*2*pi*f;                             ",
   "                                                                  ",
   "pvec=         read phase shifts from single column ASCII file     ",
   "              (radians)                                           ",
   "                                                                  ",
   " Notes:                                                           ",
   " A program that allows the user to experiment with changes in \n  ",
   " the phase spectrum of a signal.                                  ",
   "                                                                  ",
   " The original program made little real sense, e.g. the b term \n  ",
   " in the phase shift operation.   The pvec option is expecting \n  ",
   " a set of frequency dependent phase rotations."
   NULL
};

/*********************************** end self doc *********************/

/* definitions used internally */
#define PFA_MAX   720720          /* Largest allowed nfft           */

/* segy trace */
segy tr;

double pvec[SU_NFLTS];

int main( int argc, char **argv )
{
   float *rt = NULL;            /* real trace                   */

   float *amp = NULL;           /* amplitude spectra            */

   double *ph = NULL;

   register complex *ct = NULL;

   FILE *fp;

   cwp_String pvecfile = 0;

   int nt;                   /* number of points on input trace      */

   int nfft;                 /* transform length                     */

   int nf;

   float dt;                 /* sampling interval in secs            */

   float d1;                 /* output sample interval in Hz         */

   int count = 0;

   float a;

   float b;

   float c;

   float onfft;

   int i;

   double Pi = acos(-1.0);

   /* Initialize */
   initargs( argc, argv );
   requestdoc( 1 );

   /* Get info from first trace */
   if( !gettr( &tr ) )
      err( "can't get first trace" );
   nt = tr.ns;

   /* get parameters */
   /* dt is used only to set output header value d1 */
   if( !getparfloat( "dt", &dt ) ){
      dt = ( ( double ) tr.dt ) / 1000000.0;
   }

   if( !dt ) {
      dt = .004;
      warn( "dt not set, assumed to be .004" );
   }

   if( !getparfloat( "a", &a ) ){
      a = 0;
   }

   if( !getparfloat( "b", &b ) ){
      b = 1.0;
   }

   if( !getparfloat( "c", &c ) ){
      c = 0.0;
   }

   a *= Pi / 180.0;

   /* Set up pfa fft */
   nfft = npfaro( nt, 2 * nt );
   if( nfft >= SU_NFLTS || nfft >= PFA_MAX ){
      err( "Padded nt=%d--too big", nfft );
   }

   d1 = 1.0 / ( nfft * dt );
   nf = nfft / 2 + 1;
   onfft = 1.0 / nfft;

   if( getparstring( "pvec", &pvecfile ) ) {
      if( !( fp = fopen( pvecfile, "r" ) ) ) {
         perror( "unable to open pvec file" );
      } else {
         i = 0;
         while( i < nf ) {
            fscanf( fp, "%lf", &pvec[i++] );
         }
      }

      if( i != nf ) {
         sprintf( msg ,"pvec samples %d do not match data transform %d" 
                      ,i,nf );
         perror( "msg );
      }
   }
   checkpars();
   rt = ealloc1float( nfft );
   ct = ealloc1complex( nf );
   amp = ealloc1float( nf );
   ph = ealloc1double( nf );

   /* Main loop over traces */
   count = 0;
   do {
      register int i;

      /* Load trace into rt (zero-padded) */
      memcpy( ( void * ) rt, ( const void * ) &tr.data, nt * FSIZE );
      memset( ( void * ) ( rt + nt ), ( int ) '\0',
              ( nfft - nt ) * FSIZE );

      /* FFT */
      pfarc( 1, nfft, rt, ct );
      for( i = 0; i < nf; ++i ) {
         amp[i] = sqrt( ct[i].r*ct[i].r + ct[i].i*ct[i].i );
         if( amp[i] ) {
            ph[i] = atan2( ct[i].i, ct[i].r );
         } else {
            ph[i] = 0.0;
         }

         if( pvecfile ) {
            ph[i] += pvec[i];

         } else {
            ph[i] = a + b*ph[i] + c*i*d1*2*Pi;
         }
      }
      for( i = 0; i < nf; ++i ) {
         ct[i].r = amp[i] * cos( ph[i] );
         ct[i].i = amp[i] * sin( ph[i] );
      }
      pfacr( -1, nfft, ct, rt );
      for( i = 0; i < nt; ++i ){
         rt[i] *= onfft;
      }

      memcpy( ( void * ) tr.data, ( const void * ) rt, nt * FSIZE );
      puttr( &tr );

   } while( gettr( &tr ) );

   return EXIT_SUCCESS;
}
